load packages

library(tidyverse)
library(knitr)

define palettes

pal_control = c("#2E86B4", "#F2B827")
pal_outcome = wesanderson::wes_palette("Zissou1", 3, "continuous")
pal_condition = wesanderson::wes_palette("Zissou1", 6, "continuous")
pal_condition8 = wesanderson::wes_palette("Zissou1", 16, "continuous")
pal_condition8 = c(pal_condition8[1], pal_condition8[3],
                  pal_condition8[14], pal_condition8[5],
                  pal_condition8[15], pal_condition8[16],
                  pal_condition8[9], pal_condition8[12])

source data

  • Note that % fat variable was incorrect in 2012_FDES –> replaced it with the data from the EMA dataframe
source("load_data.R")

ind_diffs = ind_diffs %>%
  mutate(bmi_z = as.numeric(scale(bmi)),
         fat_z = as.numeric(scale(fat)))

standardize brain data

  • Windsorizing observations +/- 3 SD from the mean to 3 –> otherwise we lose these subs in the SCA analyses
  • Z-score within ROI and session for ROIs
  • Z-score within map, test (association or uniformity), mask (masked or unmasked) and session for PEVs
# roi distributions
betas %>%
  filter(session == "all") %>%
  select(subjectID, con, condition, control, roi, process, meanPE) %>%
  unique() %>%
  filter(!is.na(roi)) %>%
  ggplot(aes(meanPE, fill = roi)) +
  geom_density(color = NA) +
  facet_wrap(~roi, ncol = 4) +
  theme_minimal(base_size = 14) +
  theme(legend.position = "none")

# dot product distribution
dots %>%
  filter(session == "all") %>%
  select(subjectID, con, condition, control, process, test, dotProduct) %>%
  unique() %>%
  ggplot(aes(dotProduct, fill = test)) +
  geom_density(color = NA) +
  facet_wrap(~process, ncol = 3, scales = "free") +
  theme_minimal(base_size = 14) +
  theme(legend.position = c(.85, .2))

dots %>%
  filter(session == "all") %>%
  select(subjectID, con, condition, control, process, test, dotProduct, mask) %>%
  unique() %>%
  ggplot(aes(dotProduct, fill = mask, alpha = .5)) +
  geom_density(color = NA) +
  facet_wrap(~process, ncol = 3, scales = "free") +
  theme_minimal(base_size = 14) +
  theme(legend.position = c(.85, .2))

# standardize and winsorize
betas_std = betas %>%
  group_by(roi, session) %>%
  mutate(meanPE_std = scale(meanPE, center = TRUE, scale = TRUE),
         outlier = ifelse(meanPE_std > 3, "yes",
                   ifelse(meanPE_std < -3, "yes", "no")),
         meanPE_std = ifelse(meanPE_std > 3, 3,
                      ifelse(meanPE_std < -3, -3, meanPE_std))) %>%
  ungroup()

dots_std = dots %>%
  group_by(map, test, mask, session) %>%
  mutate(dotProduct_std = scale(dotProduct, center = TRUE, scale = TRUE),
         outlier = ifelse(dotProduct_std > 3, "yes",
                   ifelse(dotProduct_std < -3, "yes", "no")),
         dotProduct_std = ifelse(dotProduct_std > 3, 3,
                          ifelse(dotProduct_std < -3, -3, dotProduct_std))) %>%
  ungroup()

summarize outliers

# summarize
betas_std %>%
  filter(session == "all") %>%
  group_by(outlier) %>%
  summarize(n = n()) %>%
  spread(outlier, n) %>%
  mutate(percent = round((yes / (yes + no)) * 100, 2))
dots_std %>%
  filter(session == "all") %>%
  group_by(outlier, test, mask) %>%
  summarize(n = n()) %>%
  spread(outlier, n) %>%
  mutate(percent = round((yes / (yes + no)) * 100, 2))
# remove outlier column
betas_std = betas_std %>%
  select(-outlier)

dots_std = dots_std %>%
  select(-outlier)

merge data

dataset = full_join(betas_std, dots_std, by = c("subjectID", "con", "process", "condition", "control", "session")) %>%
  mutate(subjectID = as.character(subjectID)) %>%
  left_join(., ind_diffs, by = "subjectID") %>%
  ungroup() %>%
  mutate(subjectID = as.factor(subjectID),
         condition = as.factor(condition),
         control = as.factor(control),
         roi = as.factor(roi),
         process = as.factor(process),
         test = as.factor(test))

saveRDS(dataset, "full_dataset.RDS")
saveRDS(betas_std, "betas_dataset.RDS")
saveRDS(dots_std, "dots_dataset.RDS")

tidy data for plotting

  • Select cons generated across all sessions
  • Calculate mean across ROIs within each psychological process
  • Calculate balance score ROI = cognitive control - reward
dots_plot = dots_std %>%
  filter(mask == "unmasked") %>%
  filter(session == "all") %>%
  select(-c(map, con, mask, session, dotProduct)) %>%
  left_join(., ind_diffs) %>%
  select(-c(sample, DBIC_ID, age, gender))

betas_plot = betas_std %>%
  group_by(subjectID, process, condition, control) %>%
  mutate(meanProcessPEstd = mean(meanPE_std, na.rm = TRUE)) %>%
  select(-c(con, xyz, roi, meanPE, sdPE, meanPE_std)) %>%
  unique() %>%
  spread(process, meanProcessPEstd) %>%
  mutate(balance = cognitive_control - reward) %>%
  gather(process, meanProcessPEstd, cognitive_control, reward, value, balance) %>%
  filter(session == "all") %>%
  select(-c(session)) %>%
  left_join(., ind_diffs) %>%
  select(-c(sample, DBIC_ID, age, gender)) %>%
  ungroup()

combined_plot = dots_plot %>%
  rename("value_std" = dotProduct_std) %>%
  mutate(test = ifelse(process == "craving_regulation", "neural signature", test)) %>%
  select(subjectID, process, test, condition, control, value_std, fat, bmi) %>%
  bind_rows(betas_plot %>%
              rename("value_std" = meanProcessPEstd) %>%
              mutate(test = "ROI") %>%
              select(-contains("_z")))

tidy data for modeling

ema_tidy = ema %>%
  select(subjectID, desire_pres0.prop, desire_str.m, goal_conf.m, resist.m, enact_prop, amount_Recoded.m, hunger.m) %>%
  unique() %>%
  mutate(enact_prop_z = as.numeric(scale(enact_prop)))

model_data = dataset %>%
  filter(test == "association" & mask == "unmasked" & session == "all" & control == "nature" & condition == "food") %>%
  unique() %>%
  group_by(process, subjectID) %>%
  mutate(meanProcessPEstd = mean(meanPE_std, na.rm = TRUE),
         processPE = paste0(process, "_ROI"),
         processPEV = paste0(process, "_PEV")) %>%
  ungroup() %>%
  select(-c(xyz, meanPE, meanPE_std, sdPE, dotProduct, map, process)) %>%
  spread(processPEV, dotProduct_std) %>%
  spread(processPE, meanProcessPEstd) %>%
  group_by(subjectID) %>%
  fill(contains("PEV"), .direction = "down") %>%
  fill(contains("PEV"), .direction = "up") %>%
  fill(contains("ROI"), .direction = "down") %>%
  fill(contains("ROI"), .direction = "up") %>%
  select(-c(roi, craving_ROI, craving_regulation_ROI, test, age, gender, mask)) %>%
  unique() %>%
  mutate(balance_ROI = cognitive_control_ROI - reward_ROI)

model_data_bmi = model_data %>%
  select(-contains("fat")) %>%
  na.omit()

model_data_fat = model_data %>%
  select(-contains("bmi")) %>%
  na.omit()

model_data_ema = model_data %>%
  select(-contains("fat"), -contains("bmi")) %>%
  right_join(., ema_tidy, by = "subjectID") %>%
  na.omit()

visualize

DVs

# source raincloud plot
source("https://gist.githubusercontent.com/benmarwick/2a1bb0133ff568cbe28d/raw/fb53bd97121f7f9ce947837ef1a4c65a73bffb3f/geom_flat_violin.R")

# plot
model_data %>%
  left_join(., ema_tidy, by = "subjectID") %>%
  filter(!is.na(subjectID)) %>%
  select(subjectID, bmi, fat, enact_prop) %>%
  rename("BMI" = bmi,
         "body fat percentage" = fat,
         "enactment" = enact_prop) %>%
  gather(var, val, -subjectID) %>%
  ggplot(aes("", val, fill = var)) +
    geom_flat_violin(position = position_nudge(x = .1, y = 0), color = FALSE) +
    geom_boxplot(width = .1, outlier.shape = NA) +
    geom_point(position = position_jitter(width = .05), size = .5, alpha = .5) +
    facet_wrap(~var, scales = "free") +
    scale_fill_manual(values = pal_outcome) +
    scale_x_discrete(expand = c(.1, 0)) +
    coord_flip() +
    labs(y = "\nvalue", x = "") +
    theme_minimal(base_size = 14) +
    theme(legend.position = "none")

combined

combined_plot %>%
  filter(!process == "balance") %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature")),
         process = gsub("_", " ", process),
         test = factor(test, levels = c("neural signature", "ROI", "association", "uniformity"))) %>%
  ggplot(aes(control, value_std, fill = condition)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_wrap(~ test + process, nrow = 3) +
  scale_fill_manual(name = "", values = pal_condition) +
  labs(x = "\ncontrol condition", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

combined_plot %>%
  filter(!process == "balance") %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature")),
         process = gsub("_", " ", process),
         process = ifelse(!test == "ROI", sprintf("%s PEV", process), sprintf("%s ROI", process)),
         test = ifelse(!test %in% c("association", "uniformity"), "none", test),
         test = factor(test, levels = c("none", "association", "uniformity"))) %>%
  ggplot(aes(control, value_std, fill = process)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(test~condition) +
  scale_fill_manual(name = "", values = pal_condition8) +
  labs(x = "\ncontrol condition", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

body composition plots

ROIs

Mean parameter estimates as a function of process, stimulus type, and control condition

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(condition, meanProcessPEstd, fill = control)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(~process) +
  scale_fill_manual(name = "", values = pal_control) +
  labs(x = "\nstimulus", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(control, meanProcessPEstd, fill = condition)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(~process) +
  scale_fill_manual(name = "", values = pal_condition) +
  labs(x = "\ncontrol condition", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

Correlations between variables and BMI or % fat

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(bmi, meanProcessPEstd, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control~process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nBMI", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(fat, meanProcessPEstd, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control~process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nbody composition (% fat)", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

PEVs

Mean correspondance as a function of process, stimulus type, control condition, and test

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(condition, dotProduct_std, fill = control)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(test~process) +
  scale_fill_manual(name = "", values = pal_control) +
  labs(x = "\nstimulus", y = "mean correspondance\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(control, dotProduct_std, fill = condition)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(test~process) +
  scale_fill_manual(name = "", values = pal_condition) +
  labs(x = "\ncontrol condition", y = "mean correspondance\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

Correlations between variables and BMI or % fat

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(bmi, dotProduct_std, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control + test ~ process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nBMI", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(fat, dotProduct_std, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control + test ~ process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nbody composition (% fat)", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

EMA plots

ROIs

model_data_ema %>%
  select(subjectID, contains("ROI"), desire_pres0.prop, desire_str.m, goal_conf.m, resist.m, enact_prop, amount_Recoded.m) %>%
  unique() %>%
  gather(neuro_var, neuro_val, contains("ROI")) %>%
  mutate(neuro_var = factor(neuro_var, levels = c("balance_ROI", "cognitive_control_ROI",
                                                  "reward_ROI", "value_ROI"))) %>%
  gather(ema_var, ema_val, contains("m"), contains("prop")) %>%
  group_by(ema_var) %>%
  do({
    plot = ggplot(., aes(neuro_val, ema_val, color = neuro_var)) +
      geom_point(alpha = .5) +
      geom_smooth(method = "lm", alpha = .1, fullrange = TRUE) +
      scale_color_manual(name = "", values = pal_condition) +
      labs(x = "\n mean parameter estimate", y = paste0(.$ema_var[[1]], "\n")) +
      theme_minimal(base_size = 14)
    
    print(plot)
    data.frame()
  })

PEVs

model_data_ema %>%
  select(subjectID, contains("PEV"), desire_pres0.prop, desire_str.m, goal_conf.m, resist.m, enact_prop, amount_Recoded.m) %>%
  unique() %>%
  gather(neuro_var, neuro_val, contains("PEV")) %>%
  mutate(neuro_var = factor(neuro_var, levels = c("cognitive_control_PEV", "reward_PEV", "value_PEV",
                                                  "craving_PEV", "craving_regulation_PEV"))) %>%
  gather(ema_var, ema_val, contains("m"), contains("prop")) %>%
  group_by(ema_var) %>%
  do({
    plot = ggplot(., aes(neuro_val, ema_val, color = neuro_var)) +
      geom_point(alpha = .5) +
      geom_smooth(method = "lm", alpha = .1, fullrange = TRUE) +
      scale_color_manual(name = "", values = pal_condition[2:6]) +
      labs(x = "\n mean parameter estimate", y = paste0(.$ema_var[[1]], "\n")) +
      theme_minimal(base_size = 14)
    
    print(plot)
    data.frame()
  })

correlations among ROIs and PEVs

dots_plot %>%
  filter(!grepl("craving", process)) %>%
  left_join(., betas_plot) %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(meanProcessPEstd, dotProduct_std, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control + test ~ process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nmean standardized parameter estimate", y = "mean standardized PEV\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

PEVs v ROIs

dots_plot %>%
  unite(var, process, test, condition, control, sep = "_") %>%
  mutate(var = sprintf("dot_%s", var)) %>%
  spread(var, dotProduct_std) %>%
  left_join(., betas_plot) %>%
  unite(var, process, condition, control, sep = "_") %>%
  mutate(var = sprintf("roi_%s", var)) %>%
  spread(var, meanProcessPEstd) %>%
  select(-subjectID) %>%
  cor(., use = "pairwise.complete.obs") %>%
  reshape2::melt() %>%
  filter(!grepl("dot", Var1) & !grepl("roi", Var2)) %>%
  ggplot(aes(Var1, Var2, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c(pal_condition[1], "white", pal_condition[6]), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 1)), size = 2) +
    labs(x = "", y = "") + 
    theme(legend.text = element_text(size = 8),
          legend.position = "top",
          legend.title = element_blank(),
          axis.text = element_text(color = "black"),
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
          axis.line = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

ROIs

dots_plot %>%
  unite(var, process, test, condition, control, sep = "_") %>%
  mutate(var = sprintf("dot_%s", var)) %>%
  spread(var, dotProduct_std) %>%
  left_join(., betas_plot) %>%
  unite(var, process, condition, control, sep = "_") %>%
  mutate(var = sprintf("roi_%s", var)) %>%
  spread(var, meanProcessPEstd) %>%
  select(-subjectID) %>%
  cor(., use = "pairwise.complete.obs") %>%
  reshape2::melt() %>%
  filter(grepl("roi", Var1) & grepl("roi", Var2)) %>%
  ggplot(aes(Var1, Var2, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c(pal_condition[1], "white", pal_condition[6]), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 1)), size = 2) +
    labs(x = "", y = "") + 
    theme(legend.text = element_text(size = 8),
          legend.position = "top",
          legend.title = element_blank(),
          axis.text = element_text(color = "black"),
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
          axis.line = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

PEVs

dots_plot %>%
  unite(var, process, test, condition, control, sep = "_") %>%
  mutate(var = sprintf("dot_%s", var)) %>%
  spread(var, dotProduct_std) %>%
  left_join(., betas_plot) %>%
  unite(var, process, condition, control, sep = "_") %>%
  mutate(var = sprintf("roi_%s", var)) %>%
  spread(var, meanProcessPEstd) %>%
  select(-subjectID) %>%
  cor(., use = "pairwise.complete.obs") %>%
  reshape2::melt() %>%
  filter(grepl("dot", Var1) & grepl("dot", Var2)) %>%
  ggplot(aes(Var1, Var2, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c(pal_condition[1], "white", pal_condition[6]), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 1)), size = 2) +
    labs(x = "", y = "") + 
    theme(legend.text = element_text(size = 8),
          legend.position = "top",
          legend.title = element_blank(),
          axis.text = element_text(color = "black"),
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
          axis.line = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

run cross-validated models

descriptives

model_data %>%
  left_join(., model_data_ema) %>%
  gather(Variable, val, contains("PEV"), contains("ROI"), bmi, fat, enact_prop) %>%
  filter(!Variable == "balance_ROI") %>%
  group_by(Variable) %>%
  summarize(M = mean(val, na.rm = TRUE),
            SD = sd(val, na.rm = TRUE)) %>%
  filter(!is.na(Variable)) %>%
  mutate(Variable = gsub("_", " ", Variable),
         Variable = ifelse(Variable == "bmi", "BMI",
                    ifelse(Variable == "fat", "Body fat percentage",
                    ifelse(Variable == "enact prop", "Enactment", Hmisc::capitalize(Variable)))),
         Variable = factor(Variable, levels = c("BMI", "Body fat percentage", "Enactment",
                                                "Reward ROI", "Reward PEV",
                                                "Value ROI", "Value PEV",
                                                "Craving regulation PEV", "Craving PEV",
                                                "Cognitive control ROI", "Cognitive control PEV"))) %>%
  arrange(Variable) %>%
  kable(format = "pandoc", digits = 2)
Variable M SD
BMI 23.19 3.36
Body fat percentage 25.43 8.00
Enactment 0.47 0.22
Reward ROI 0.06 0.50
Reward PEV 0.07 0.77
Value ROI 0.40 0.68
Value PEV 0.28 0.77
Craving regulation PEV 0.40 0.67
Craving PEV -0.13 0.76
Cognitive control ROI -0.13 0.50
Cognitive control PEV -0.11 0.78
cor_mat = model_data %>%
  left_join(., model_data_ema) %>%
  select(bmi, fat, enact_prop, contains("ROI"), contains("PEV"), -balance_ROI) %>%
  gather(var, val, contains("PEV"), contains("ROI"), bmi, fat, enact_prop) %>%
  mutate(var = gsub("_", " ", var),
         var = ifelse(var == "bmi", "BMI",
               ifelse(var == "fat", "Body fat percentage",
               ifelse(var == "enact prop", "Enactment", Hmisc::capitalize(var)))),
         var = factor(var, levels = c("BMI", "Body fat percentage", "Enactment",
                                      "Reward ROI", "Reward PEV",
                                      "Value ROI", "Value PEV",
                                      "Craving regulation PEV", "Craving PEV",
                                      "Cognitive control ROI", "Cognitive control PEV"))) %>%
  arrange(var) %>%
  spread(var, val) %>%
  ungroup() %>%
  select(-subjectID) %>%
  as.matrix()

cors = Hmisc::rcorr(cor_mat)

cors$r %>%
  as.data.frame() %>%
  kable(format = "pandoc", digits = 2)
BMI Body fat percentage Enactment Reward ROI Reward PEV Value ROI Value PEV Craving regulation PEV Craving PEV Cognitive control ROI Cognitive control PEV
BMI 1.00 0.66 -0.13 -0.12 -0.02 0.06 0.00 0.10 -0.02 -0.02 -0.06
Body fat percentage 0.66 1.00 -0.20 -0.07 -0.07 0.01 -0.04 0.02 -0.07 -0.03 -0.06
Enactment -0.13 -0.20 1.00 0.21 0.31 0.23 0.29 0.03 0.21 0.01 0.05
Reward ROI -0.12 -0.07 0.21 1.00 0.62 0.30 0.60 0.08 0.47 0.31 0.41
Reward PEV -0.02 -0.07 0.31 0.62 1.00 0.58 0.91 0.19 0.81 0.47 0.59
Value ROI 0.06 0.01 0.23 0.30 0.58 1.00 0.65 0.28 0.51 0.06 0.16
Value PEV 0.00 -0.04 0.29 0.60 0.91 0.65 1.00 0.15 0.76 0.30 0.40
Craving regulation PEV 0.10 0.02 0.03 0.08 0.19 0.28 0.15 1.00 0.24 0.21 0.23
Craving PEV -0.02 -0.07 0.21 0.47 0.81 0.51 0.76 0.24 1.00 0.45 0.54
Cognitive control ROI -0.02 -0.03 0.01 0.31 0.47 0.06 0.30 0.21 0.45 1.00 0.91
Cognitive control PEV -0.06 -0.06 0.05 0.41 0.59 0.16 0.40 0.23 0.54 0.91 1.00
cors$P %>%
  as.data.frame() %>%
  mutate_if(is.numeric, funs(ifelse(between(., .05, .1), "†",
                             ifelse(between(., .01, .05), "*", 
                             ifelse(between(., .001, .01), "**",
                             ifelse(. < .001, "***", "")))))) %>%
  kable(format = "pandoc", digits = 2)
BMI Body fat percentage Enactment Reward ROI Reward PEV Value ROI Value PEV Craving regulation PEV Craving PEV Cognitive control ROI Cognitive control PEV
NA ***
*** NA
NA * ** * ** *
* NA *** *** *** *** *** ***
** *** NA *** *** ** *** *** ***
* *** *** NA *** *** *** **
** *** *** *** NA * *** *** ***
** *** * NA *** *** ***
* *** *** *** *** *** NA *** ***
*** *** *** *** *** NA ***
*** *** ** *** *** *** *** NA

t-tests

model_data %>%
  select(-contains("z"), -balance_ROI) %>%
  left_join(., select(model_data_ema, subjectID, enact_prop)) %>%
  ungroup() %>%
  select_if(is.numeric) %>%
  map_df(~ broom::tidy(t.test(.x, mu = 0, na.action = "na.omit")), .id = "x") %>%
  rename("Variable" = x,
         "t" = statistic,
         "df" = parameter,
         "p" = p.value) %>%
  mutate(`Mdiff [95% CI]` = sprintf("%.02f [%.02f, %.02f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         Variable = gsub("enact_prop", "Enactment", Variable),
         Variable = gsub("fat", "Body fat percentage", Variable),
         Variable = gsub("bmi", "BMI", Variable),
         Variable = gsub("_", " ", Variable),
         Variable = factor(Variable, levels = c("BMI", "Body fat percentage", "Enactment",
                                                "reward ROI", "reward PEV",
                                                "value ROI", "value PEV",
                                                "craving regulation PEV", "craving PEV",
                                                "cognitive control ROI", "cognitive control PEV"))) %>%
  select(Variable, `Mdiff [95% CI]`, t, df, p) %>%
  arrange(Variable) %>%
  kable(format = "pandoc", digits = 2)
Variable Mdiff [95% CI] t df p
BMI 23.19 [22.77, 23.60] 109.65 252 < .001
Body fat percentage 25.43 [24.44, 26.43] 50.28 249 < .001
Enactment 0.47 [0.42, 0.51] 20.58 96 < .001
reward ROI 0.06 [-0.00, 0.12] 1.90 261 .059
reward PEV 0.07 [-0.02, 0.17] 1.51 261 .132
value ROI 0.40 [0.31, 0.48] 9.38 261 < .001
value PEV 0.28 [0.18, 0.37] 5.84 261 < .001
craving regulation PEV 0.40 [0.31, 0.48] 9.52 261 < .001
craving PEV -0.13 [-0.22, -0.04] -2.74 261 .007
cognitive control ROI -0.13 [-0.19, -0.07] -4.14 261 < .001
cognitive control PEV -0.11 [-0.20, -0.01] -2.20 261 .029

specify model variables

options(na.action = "na.fail")
data_ctrl = caret::trainControl(method = "repeatedcv", number = 5, repeats = 5)

BMI

ROI models

bmi_roi_full = lm(bmi_z ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data_bmi)
(bmi_roi_mods = MuMIn::dredge(bmi_roi_full))
MuMIn::get.models(bmi_roi_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()
best_bmi_roi = caret::train(bmi_z ~ reward_ROI + value_ROI,
                     data = model_data_bmi,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

PEV models

bmi_pev_full = lm(bmi_z ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data_bmi)
(bmi_pev_mods = MuMIn::dredge(bmi_pev_full))
MuMIn::get.models(bmi_pev_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()
best_bmi_pev = caret::train(bmi_z ~ craving_regulation_PEV,
                     data = model_data_bmi,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

combined

bmi_null = lm(bmi_z ~ 1, data = model_data_bmi)
summary(bmi_null)
## 
## Call:
## lm(formula = bmi_z ~ 1, data = model_data_bmi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2100 -0.6542 -0.0558  0.4827  4.6115 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.01306    0.06326   0.206    0.837
## 
## Residual standard error: 1.006 on 252 degrees of freedom
bmi_combined = caret::train(bmi_z ~ reward_ROI + value_ROI + craving_regulation_PEV,
                     data = model_data_bmi,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

(bmi_aic = AIC(bmi_null, best_bmi_roi$finalModel, best_bmi_pev$finalModel, bmi_combined$finalModel) %>%
    rownames_to_column() %>%
    extract(rowname, "model", ".*bmi_(.*)\\$.*") %>%
    mutate(model = toupper(model),
           model = ifelse(is.na(model), "Null", model)) %>%
  arrange(AIC))

summarize the overall best fitting

# ROI
summary(best_bmi_roi$finalModel)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1355 -0.6606 -0.0899  0.4662  4.5052 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.03113    0.07236  -0.430   0.6674  
## reward_ROI  -0.30523    0.12990  -2.350   0.0196 *
## value_ROI    0.15583    0.09654   1.614   0.1077  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9976 on 250 degrees of freedom
## Multiple R-squared:  0.02502,    Adjusted R-squared:  0.01722 
## F-statistic: 3.207 on 2 and 250 DF,  p-value: 0.04214
best_bmi_roi
## Linear Regression 
## 
## 253 samples
##   2 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 202, 201, 203, 203, 203, 203, ... 
## Resampling results:
## 
##   RMSE       Rsquared    MAE      
##   0.9949716  0.03211436  0.7469028
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE
sd(best_bmi_roi$resample$Rsquared)
## [1] 0.03262441
sd(best_bmi_roi$resample$RMSE)
## [1] 0.1099513
# combined
summary(bmi_combined$finalModel)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1994 -0.6654 -0.0886  0.5098  4.4311 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)  
## (Intercept)            -0.07322    0.07768  -0.942   0.3469  
## reward_ROI             -0.30708    0.12961  -2.369   0.0186 *
## value_ROI               0.11695    0.09989   1.171   0.2428  
## craving_regulation_PEV  0.14370    0.09792   1.467   0.1435  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9953 on 249 degrees of freedom
## Multiple R-squared:  0.03338,    Adjusted R-squared:  0.02173 
## F-statistic: 2.866 on 3 and 249 DF,  p-value: 0.03725
bmi_combined
## Linear Regression 
## 
## 253 samples
##   3 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 200, 203, 204, 203, 202, 204, ... 
## Resampling results:
## 
##   RMSE       Rsquared    MAE      
##   0.9955154  0.03419623  0.7523059
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE
sd(bmi_combined$resample$Rsquared)
## [1] 0.03657507
sd(bmi_combined$resample$RMSE)
## [1] 0.1027421

% fat

ROI models

fat_roi_full = lm(fat_z ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data_fat)
(fat_roi_mods = MuMIn::dredge(fat_roi_full))
MuMIn::get.models(fat_roi_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()
best_fat_roi = caret::train(fat_z ~ reward_ROI,
                     data = model_data_fat,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

PEV models

fat_pev_full = lm(fat_z ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data_fat)
(fat_pev_mods = MuMIn::dredge(fat_pev_full))
MuMIn::get.models(fat_pev_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()
best_fat_pev = caret::train(fat_z ~ craving_PEV,
                     data = model_data_fat,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

combined

fat_null = lm(fat_z ~ 1, data = model_data_fat)

fat_combined = caret::train(fat_z ~ reward_ROI + craving_PEV,
                     data = model_data_fat,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

(fat_aic = AIC(fat_null, best_fat_roi$finalModel, best_fat_pev$finalModel, fat_combined$finalModel) %>%
    rownames_to_column() %>%
    extract(rowname, "model", ".*fat_(.*)\\$.*") %>%
    mutate(model = toupper(model),
           model = ifelse(is.na(model), "Null", model)) %>%
    arrange(AIC))

summarize the overall best fitting

summary(fat_null)
## 
## Call:
## lm(formula = fat_z ~ 1, data = model_data_fat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9531 -0.5945  0.1396  0.6457  2.4201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.01679    0.06321   0.266    0.791
## 
## Residual standard error: 0.9994 on 249 degrees of freedom

EMA enactment proportion

ROI models

enact_prop_roi_full = lm(enact_prop_z ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data_ema)
(enact_prop_roi_mods = MuMIn::dredge(enact_prop_roi_full))
MuMIn::get.models(enact_prop_roi_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()
best_enact_prop_roi = caret::train(enact_prop_z ~ value_ROI,
                     data = model_data_ema,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

PEV models

enact_prop_pev_full = lm(enact_prop_z ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data_ema)
(enact_prop_pev_mods = MuMIn::dredge(enact_prop_pev_full))
MuMIn::get.models(enact_prop_pev_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()
best_enact_prop_pev = caret::train(enact_prop_z ~ reward_PEV,
                     data = model_data_ema,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

combined

enact_prop_null = lm(enact_prop_z ~ 1, data = model_data_ema)

enact_prop_combined = caret::train(enact_prop_z ~ value_ROI + reward_PEV,
                     data = model_data_ema,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

(enact_aic = AIC(enact_prop_null, best_enact_prop_roi$finalModel, best_enact_prop_pev$finalModel, enact_prop_combined$finalModel) %>%
    rownames_to_column() %>%
    extract(rowname, "model", ".*prop_(.*)\\$.*") %>%
    mutate(model = toupper(model),
           model = ifelse(is.na(model), "Null", model)) %>%
  arrange(AIC))

summarize the overall best fitting

summary(best_enact_prop_pev$finalModel)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.17754 -0.58154 -0.01067  0.53885  2.40418 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.09627    0.09985  -0.964  0.33745   
## reward_PEV   0.37158    0.11641   3.192  0.00192 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9697 on 95 degrees of freedom
## Multiple R-squared:  0.09686,    Adjusted R-squared:  0.08735 
## F-statistic: 10.19 on 1 and 95 DF,  p-value: 0.001917
best_enact_prop_pev
## Linear Regression 
## 
## 97 samples
##  1 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 77, 79, 77, 78, 77, 78, ... 
## Resampling results:
## 
##   RMSE       Rsquared   MAE      
##   0.9722222  0.1498315  0.7677445
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE
sd(best_enact_prop_pev$resample$Rsquared)
## [1] 0.1591469
sd(best_enact_prop_pev$resample$RMSE)
## [1] 0.1075649

tables

bmi

bmi_fit = data.frame(model = c("ROI", "PEV", "COMBINED"),
                     r2 = c(mean(best_bmi_roi$resample$Rsquared), 
                            mean(best_bmi_pev$resample$Rsquared), 
                            mean(bmi_combined$resample$Rsquared)),
                     r2_sd = c(sd(best_bmi_roi$resample$Rsquared), 
                               sd(best_bmi_pev$resample$Rsquared), 
                               sd(bmi_combined$resample$Rsquared)),
                     RMSE = c(mean(best_bmi_roi$resample$RMSE), 
                              mean(best_bmi_pev$resample$RMSE), 
                              mean(bmi_combined$resample$RMSE)),
                     RMSE_sd = c(sd(best_bmi_roi$resample$RMSE), 
                                 sd(best_bmi_pev$resample$RMSE), 
                                 sd(bmi_combined$resample$RMSE)))

bmi_table = broom::tidy(best_bmi_roi$finalModel, conf.int = TRUE) %>% 
  mutate(model = "ROI", 
         mod_sig = ifelse(pf(summary(best_bmi_roi$finalModel)$fstatistic[1], 
                             summary(best_bmi_roi$finalModel)$fstatistic[2], 
                             summary(best_bmi_roi$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", "")) %>%
  bind_rows(broom::tidy(best_bmi_pev$finalModel, conf.int = TRUE) %>% 
              mutate(model = "PEV",
                     mod_sig = ifelse(pf(summary(best_bmi_pev$finalModel)$fstatistic[1], 
                             summary(best_bmi_pev$finalModel)$fstatistic[2], 
                             summary(best_bmi_pev$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(bmi_combined$finalModel, conf.int = TRUE) %>% 
              mutate(model = "COMBINED",
                     mod_sig = ifelse(pf(summary(bmi_combined$finalModel)$fstatistic[1], 
                             summary(bmi_combined$finalModel)$fstatistic[2], 
                             summary(bmi_combined$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(bmi_null, conf.int = TRUE) %>% mutate(model = "Null", mod_sig = "")) %>%
  rename("SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept", term),
         term = gsub("_", " ", term),
         term = Hmisc::capitalize(term)) %>%
  left_join(., bmi_aic) %>%
  left_join(., bmi_fit) %>%
  mutate(`r2 (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", r2, r2_sd), "--"),
         `RMSE (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", RMSE, RMSE_sd), "--")) %>%
  unite(model, model, mod_sig, sep = "") %>%
  select(model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p) %>%
  arrange(AIC) %>%
  mutate(model = gsub("COMBINED", "Combined", model),
         outcome = "BMI") %>%
    select(outcome, model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p)

bmi_table %>%
  kable(format = "pandoc", digits = 2)
outcome model AIC r2 (SD) RMSE (SD) term b [95% CI] SE t p
BMI Combined* 721.56 0.03 (0.04) 1.00 (0.10) Intercept -0.07 [-0.23, 0.08] 0.08 -0.94 .347
BMI Combined* 721.56 0.03 (0.04) 1.00 (0.10) Reward ROI -0.31 [-0.56, -0.05] 0.13 -2.37 .019
BMI Combined* 721.56 0.03 (0.04) 1.00 (0.10) Value ROI 0.12 [-0.08, 0.31] 0.10 1.17 .243
BMI Combined* 721.56 0.03 (0.04) 1.00 (0.10) Craving regulation PEV 0.14 [-0.05, 0.34] 0.10 1.47 .144
BMI ROI* 721.74 0.03 (0.03) 0.99 (0.11) Intercept -0.03 [-0.17, 0.11] 0.07 -0.43 .667
BMI ROI* 721.74 0.03 (0.03) 0.99 (0.11) Reward ROI -0.31 [-0.56, -0.05] 0.13 -2.35 .020
BMI ROI* 721.74 0.03 (0.03) 0.99 (0.11) Value ROI 0.16 [-0.03, 0.35] 0.10 1.61 .108
BMI PEV 723.45 0.03 (0.03) 1.00 (0.13) Intercept -0.05 [-0.19, 0.10] 0.07 -0.67 .507
BMI PEV 723.45 0.03 (0.03) 1.00 (0.13) Craving regulation PEV 0.16 [-0.03, 0.34] 0.09 1.64 .102
BMI Null 724.15 Intercept 0.01 [-0.11, 0.14] 0.06 0.21 .837

% fat

fat_fit = data.frame(model = c("ROI", "PEV", "COMBINED"),
                     r2 = c(mean(best_fat_roi$resample$Rsquared), 
                            mean(best_fat_pev$resample$Rsquared), 
                            mean(fat_combined$resample$Rsquared)),
                     r2_sd = c(sd(best_fat_roi$resample$Rsquared), 
                               sd(best_fat_pev$resample$Rsquared), 
                               sd(fat_combined$resample$Rsquared)),
                     RMSE = c(mean(best_fat_roi$resample$RMSE), 
                              mean(best_fat_pev$resample$RMSE), 
                              mean(fat_combined$resample$RMSE)),
                     RMSE_sd = c(sd(best_fat_roi$resample$RMSE), 
                                 sd(best_fat_pev$resample$RMSE), 
                                 sd(fat_combined$resample$RMSE)))

fat_table = broom::tidy(best_fat_roi$finalModel, conf.int = TRUE) %>% 
  mutate(model = "ROI", 
         mod_sig = ifelse(pf(summary(best_fat_roi$finalModel)$fstatistic[1], 
                             summary(best_fat_roi$finalModel)$fstatistic[2], 
                             summary(best_fat_roi$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", "")) %>%
  bind_rows(broom::tidy(best_fat_pev$finalModel, conf.int = TRUE) %>% 
              mutate(model = "PEV",
                     mod_sig = ifelse(pf(summary(best_fat_pev$finalModel)$fstatistic[1], 
                             summary(best_fat_pev$finalModel)$fstatistic[2], 
                             summary(best_fat_pev$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(fat_combined$finalModel, conf.int = TRUE) %>% 
              mutate(model = "COMBINED",
                     mod_sig = ifelse(pf(summary(fat_combined$finalModel)$fstatistic[1], 
                             summary(fat_combined$finalModel)$fstatistic[2], 
                             summary(fat_combined$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(fat_null, conf.int = TRUE) %>% mutate(model = "Null", mod_sig = "")) %>%
  rename("SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept", term),
         term = gsub("_", " ", term),
         term = Hmisc::capitalize(term)) %>%
  left_join(., fat_aic) %>%
  left_join(., fat_fit) %>%
  mutate(`r2 (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", r2, r2_sd), "--"),
         `RMSE (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", RMSE, RMSE_sd), "--")) %>%
  unite(model, model, mod_sig, sep = "") %>%
  select(model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p) %>%
  arrange(AIC) %>%
  mutate(model = gsub("COMBINED", "Combined", model),
         outcome = "Body fat") %>%
    select(outcome, model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p)

fat_table %>%
  kable(format = "pandoc", digits = 2)
outcome model AIC r2 (SD) RMSE (SD) term b [95% CI] SE t p
Body fat Null 712.18 Intercept 0.02 [-0.11, 0.14] 0.06 0.27 .791
Body fat ROI 712.92 0.02 (0.02) 1.00 (0.07) Intercept 0.02 [-0.10, 0.15] 0.06 0.38 .706
Body fat ROI 712.92 0.02 (0.02) 1.00 (0.07) Reward ROI -0.14 [-0.39, 0.11] 0.13 -1.12 .265
Body fat PEV 712.94 0.03 (0.03) 1.00 (0.08) Intercept 0.00 [-0.12, 0.13] 0.06 0.05 .960
Body fat PEV 712.94 0.03 (0.03) 1.00 (0.08) Craving PEV -0.09 [-0.26, 0.07] 0.08 -1.11 .268
Body fat Combined 714.49 0.02 (0.03) 1.01 (0.06) Intercept 0.01 [-0.12, 0.14] 0.07 0.19 .849
Body fat Combined 714.49 0.02 (0.03) 1.01 (0.06) Reward ROI -0.09 [-0.38, 0.19] 0.14 -0.66 .508
Body fat Combined 714.49 0.02 (0.03) 1.01 (0.06) Craving PEV -0.06 [-0.25, 0.13] 0.10 -0.65 .516

EMA enactment

enact_prop_fit = data.frame(model = c("ROI", "PEV", "COMBINED"),
                     r2 = c(mean(best_enact_prop_roi$resample$Rsquared), 
                            mean(best_enact_prop_pev$resample$Rsquared), 
                            mean(enact_prop_combined$resample$Rsquared)),
                     r2_sd = c(sd(best_enact_prop_roi$resample$Rsquared), 
                               sd(best_enact_prop_pev$resample$Rsquared), 
                               sd(enact_prop_combined$resample$Rsquared)),
                     RMSE = c(mean(best_enact_prop_roi$resample$RMSE), 
                              mean(best_enact_prop_pev$resample$RMSE), 
                              mean(enact_prop_combined$resample$RMSE)),
                     RMSE_sd = c(sd(best_enact_prop_roi$resample$RMSE), 
                                 sd(best_enact_prop_pev$resample$RMSE), 
                                 sd(enact_prop_combined$resample$RMSE)))

enact_table = broom::tidy(best_enact_prop_roi$finalModel, conf.int = TRUE) %>% 
  mutate(model = "ROI", 
         mod_sig = ifelse(pf(summary(best_enact_prop_roi$finalModel)$fstatistic[1], 
                             summary(best_enact_prop_roi$finalModel)$fstatistic[2], 
                             summary(best_enact_prop_roi$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", "")) %>%
  bind_rows(broom::tidy(best_enact_prop_pev$finalModel, conf.int = TRUE) %>% 
              mutate(model = "PEV",
                     mod_sig = ifelse(pf(summary(best_enact_prop_pev$finalModel)$fstatistic[1], 
                             summary(best_enact_prop_pev$finalModel)$fstatistic[2], 
                             summary(best_enact_prop_pev$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(enact_prop_combined$finalModel, conf.int = TRUE) %>% 
              mutate(model = "COMBINED",
                     mod_sig = ifelse(pf(summary(enact_prop_combined$finalModel)$fstatistic[1], 
                             summary(enact_prop_combined$finalModel)$fstatistic[2], 
                             summary(enact_prop_combined$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(enact_prop_null, conf.int = TRUE) %>% mutate(model = "Null", mod_sig = "")) %>%
  rename("SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept", term),
         term = gsub("_", " ", term),
         term = Hmisc::capitalize(term)) %>%
  left_join(., enact_aic) %>%
  left_join(., enact_prop_fit) %>%
  mutate(`r2 (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", r2, r2_sd), "--"),
         `RMSE (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", RMSE, RMSE_sd), "--")) %>%
  unite(model, model, mod_sig, sep = "") %>%
  select(model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p) %>%
  arrange(AIC) %>%
  mutate(model = gsub("COMBINED", "Combined", model),
         outcome = "Enactment") %>%
    select(outcome, model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p)

enact_table %>%
  kable(format = "pandoc", digits = 2)
outcome model AIC r2 (SD) RMSE (SD) term b [95% CI] SE t p
Enactment PEV* 273.28 0.15 (0.16) 0.97 (0.11) Intercept -0.10 [-0.29, 0.10] 0.10 -0.96 .337
Enactment PEV* 273.28 0.15 (0.16) 0.97 (0.11) Reward PEV 0.37 [0.14, 0.60] 0.12 3.19 .002
Enactment Combined* 274.60 0.12 (0.13) 0.97 (0.13) Intercept -0.15 [-0.38, 0.09] 0.12 -1.25 .215
Enactment Combined* 274.60 0.12 (0.13) 0.97 (0.13) Value ROI 0.13 [-0.19, 0.46] 0.16 0.81 .419
Enactment Combined* 274.60 0.12 (0.13) 0.97 (0.13) Reward PEV 0.32 [0.05, 0.58] 0.14 2.33 .022
Enactment ROI* 278.06 0.10 (0.09) 0.99 (0.14) Intercept -0.19 [-0.43, 0.05] 0.12 -1.59 .115
Enactment ROI* 278.06 0.10 (0.09) 0.99 (0.14) Value ROI 0.33 [0.04, 0.61] 0.14 2.27 .026
Enactment Null 281.16 Intercept -0.04 [-0.25, 0.16] 0.10 -0.42 .676

combined

bind_rows(bmi_table, fat_table, enact_table) %>%
  kable(format = "pandoc", digits = 2)
outcome model AIC r2 (SD) RMSE (SD) term b [95% CI] SE t p
BMI Combined* 721.56 0.03 (0.04) 1.00 (0.10) Intercept -0.07 [-0.23, 0.08] 0.08 -0.94 .347
BMI Combined* 721.56 0.03 (0.04) 1.00 (0.10) Reward ROI -0.31 [-0.56, -0.05] 0.13 -2.37 .019
BMI Combined* 721.56 0.03 (0.04) 1.00 (0.10) Value ROI 0.12 [-0.08, 0.31] 0.10 1.17 .243
BMI Combined* 721.56 0.03 (0.04) 1.00 (0.10) Craving regulation PEV 0.14 [-0.05, 0.34] 0.10 1.47 .144
BMI ROI* 721.74 0.03 (0.03) 0.99 (0.11) Intercept -0.03 [-0.17, 0.11] 0.07 -0.43 .667
BMI ROI* 721.74 0.03 (0.03) 0.99 (0.11) Reward ROI -0.31 [-0.56, -0.05] 0.13 -2.35 .020
BMI ROI* 721.74 0.03 (0.03) 0.99 (0.11) Value ROI 0.16 [-0.03, 0.35] 0.10 1.61 .108
BMI PEV 723.45 0.03 (0.03) 1.00 (0.13) Intercept -0.05 [-0.19, 0.10] 0.07 -0.67 .507
BMI PEV 723.45 0.03 (0.03) 1.00 (0.13) Craving regulation PEV 0.16 [-0.03, 0.34] 0.09 1.64 .102
BMI Null 724.15 Intercept 0.01 [-0.11, 0.14] 0.06 0.21 .837
Body fat Null 712.18 Intercept 0.02 [-0.11, 0.14] 0.06 0.27 .791
Body fat ROI 712.92 0.02 (0.02) 1.00 (0.07) Intercept 0.02 [-0.10, 0.15] 0.06 0.38 .706
Body fat ROI 712.92 0.02 (0.02) 1.00 (0.07) Reward ROI -0.14 [-0.39, 0.11] 0.13 -1.12 .265
Body fat PEV 712.94 0.03 (0.03) 1.00 (0.08) Intercept 0.00 [-0.12, 0.13] 0.06 0.05 .960
Body fat PEV 712.94 0.03 (0.03) 1.00 (0.08) Craving PEV -0.09 [-0.26, 0.07] 0.08 -1.11 .268
Body fat Combined 714.49 0.02 (0.03) 1.01 (0.06) Intercept 0.01 [-0.12, 0.14] 0.07 0.19 .849
Body fat Combined 714.49 0.02 (0.03) 1.01 (0.06) Reward ROI -0.09 [-0.38, 0.19] 0.14 -0.66 .508
Body fat Combined 714.49 0.02 (0.03) 1.01 (0.06) Craving PEV -0.06 [-0.25, 0.13] 0.10 -0.65 .516
Enactment PEV* 273.28 0.15 (0.16) 0.97 (0.11) Intercept -0.10 [-0.29, 0.10] 0.10 -0.96 .337
Enactment PEV* 273.28 0.15 (0.16) 0.97 (0.11) Reward PEV 0.37 [0.14, 0.60] 0.12 3.19 .002
Enactment Combined* 274.60 0.12 (0.13) 0.97 (0.13) Intercept -0.15 [-0.38, 0.09] 0.12 -1.25 .215
Enactment Combined* 274.60 0.12 (0.13) 0.97 (0.13) Value ROI 0.13 [-0.19, 0.46] 0.16 0.81 .419
Enactment Combined* 274.60 0.12 (0.13) 0.97 (0.13) Reward PEV 0.32 [0.05, 0.58] 0.14 2.33 .022
Enactment ROI* 278.06 0.10 (0.09) 0.99 (0.14) Intercept -0.19 [-0.43, 0.05] 0.12 -1.59 .115
Enactment ROI* 278.06 0.10 (0.09) 0.99 (0.14) Value ROI 0.33 [0.04, 0.61] 0.14 2.27 .026
Enactment Null 281.16 Intercept -0.04 [-0.25, 0.16] 0.10 -0.42 .676

VIF

# BMI
car::vif(bmi_combined$finalModel)
##             reward_ROI              value_ROI craving_regulation_PEV 
##               1.101738               1.185011               1.085275
# %fat
car::vif(fat_combined$finalModel)
##  reward_ROI craving_PEV 
##    1.303633    1.303633
# enactment
car::vif(enact_prop_combined$finalModel)
##  value_ROI reward_PEV 
##   1.348297   1.348297